!>\brief
!! This module collects the computation of the local matrices for the
!! LDG method.
!!
!! \n
!!
!! This module is similar to mod_ldgh_locmat, since the local matrices
!! of the hybridized and not hybridized local discontinuous Galerkin
!! formulations are similar. The main difference is that in this
!! module we consider both Dirichlet and Neumann boundary conditions,
!! since they enforced weekly, while in the LDG-H method Dirichlet
!! boundary conditions are enforced strongly on the Lagrange
!! multiplier.
!!
!! \note There are no consistency checks on the coefficients returned
!! by the \c coeff* functions of \c mod_testcases, which thus have to
!! be correct. In particular, the diffusion tensor must be symmetric.
!<----------------------------------------------------------------------
module mod_ldg_locmat

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_linal, only: &
   mod_linal_initialized, &
   invmat_chol
 
 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_e, t_s,        &
   affmap, diameter

 use mod_bcs, only: &
   mod_bcs_initialized, &
   t_b_v,   t_b_s,   t_b_e,   &
   p_t_b_v, p_t_b_s, p_t_b_e, &
   b_dir, b_neu
 
 use mod_tau, only: &
   mod_tau_initialized, &
   t_eltau

 use mod_testcases, only: &
   mod_testcases_initialized, &
   coeff_diff,  &
   coeff_adv,   &
   coeff_re,    &
   coeff_f,     &
   coeff_dir,   &
   coeff_neu,   &
   coeff_rob,   &
   coeff_alpha

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_ldg_locmat_constructor, &
   mod_ldg_locmat_destructor,  &
   mod_ldg_locmat_initialized, &
   ldg_locmat_e, ldg_locmat_s, &
   ldg_locmat_bs, t_bcs_error

 private

!-----------------------------------------------------------------------

! Module types and parameters
 
 !> This type is used to indicate that an error has occurred in the
 !! computation of the boundary condition
 !<
 type t_bcs_error
   logical :: lerr = .false.
   character(len=1000) :: message
 end type t_bcs_error

! Module variables

 ! public members
 logical, protected ::               &
   mod_ldg_locmat_initialized = .false.
 ! private members
 real(wp), allocatable :: cckr(:,:), ccsr(:,:), ddsr(:,:)
 character(len=*), parameter :: &
   this_mod_name = 'mod_ldg_locmat'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_ldg_locmat_constructor(base)
  type(t_base), intent(in) :: base

  integer :: l, i, j
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_linal_initialized.eqv..false.)    .or. &
       (mod_master_el_initialized.eqv..false.).or. &
       (mod_base_initialized.eqv..false.)     .or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_bcs_initialized.eqv..false.)      .or. &
       (mod_tau_initialized.eqv..false.)      .or. &
       (mod_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_ldg_locmat_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   !> Precompute some local matrices.

   allocate( cckr(base%pk,base%mk) )
   cckr = 0.0_wp
   do l=1,base%m
     do i=1,base%pk
       do j=1,base%mk
         cckr(i,j) = cckr(i,j) + base%wg(l) * &
           base%divo(j,l)*base%p(i,l)
       enddo
     enddo
   enddo

   allocate( ccsr(2*base%pk,2*base%mk) )
   ! Implementation missing

   allocate( ddsr(2*base%pk,2*base%pk) )
   ! Implementation missing

   mod_ldg_locmat_initialized = .true.
 end subroutine mod_ldg_locmat_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_ldg_locmat_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_ldg_locmat_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(cckr,ccsr,ddsr)

   mod_ldg_locmat_initialized = .false.
 end subroutine mod_ldg_locmat_destructor

!-----------------------------------------------------------------------
 
 !> Compute the local matrices and right hand side (element
 !! contributions).
 !!
 !! The local matrices correspond to the following terms (there are
 !! no advective fluxes yet, however the corresponding matrces are
 !! computed and set to zero for simplicity):
 !! \f{displaymath}{
 !!  A^K_{ij} = \int_K \omega_j\cdot \mu^{-1} \omega_i\,dx, \qquad
 !!  B^K_{ij} = -C^K_{ji},
 !! \f}
 !! \f{displaymath}{
 !!  C^K_{ij} = \int_K \nabla\cdot\omega_j\,\varphi_i \,dx,
 !!   \qquad
 !!  D^K_{ij} = 0,
 !!   \qquad
 !!  f^{K,2}_{i} = \int_K f\varphi_i\,dx.
 !! \f}
 !! Notice that boundary condition terms are only considered in side
 !! local matrix computations.
 !<
 pure subroutine ldg_locmat_e(aak , bbk ,      &
                              cck , ddk , fk , &
                              base,e)
  real(wp), intent(out) :: &
    aak(:,:) , bbk(:,:) ,       &
    cck(:,:) , ddk(:,:) , fk(:)
  type(t_base), intent(in) :: base
  type(t_e), intent(in) :: e

  integer :: l, i, j
  real(wp) :: &
    mu(e%m,e%m,base%m), mui(e%m,e%m,base%m), f(base%m), &
    mui_b(e%m,e%d), bt_mui_b(e%d,e%d), mu_o(e%d), swg,  &
    xg(e%m,base%m)
  character(len=*), parameter :: &
    this_sub_name = 'ldg_locmat_e'

   ! Gaussian nodes in physical space
   xg = affmap(e,base%xig)

   ! problem coefficients
   mu    = coeff_diff(xg)
   do l=1,base%m
     call invmat_chol( mu(:,:,l) , mui(:,:,l) )
   enddo
   f     = coeff_f(xg)

   !--------------------------------------------------------------------
   ! First equation: compute aak

   ! first the volume terms
   aak = 0.0_wp
   do l=1,base%m
     ! scaled quad weight
     swg = base%wg(l)/e%det_b
     ! mui_b = mui * b
     mui_b = matmul( mui(:,:,l) , e%b )
     ! bt_mui_b = b^T * mui * b
     bt_mui_b = matmul( transpose(e%b) , mui_b )
     do i=1,base%mk
       ! mu_o = b^T * mui * b * o
       mu_o = matmul( bt_mui_b , base%o(:,i,l) )
       do j=1,base%mk
         aak(i,j) = aak(i,j) +                     &
           swg * dot_product( base%o(:,j,l) , mu_o )
       enddo
     enddo
   enddo
   ! since we know that aak is symmetric, we can fix roundoff errors
   aak = 0.5_wp * (aak+transpose(aak))
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Second equation: compute cck and ddk
   cck = cckr
   ddk = 0.0_wp
   fk  = 0.0_wp
   do l=1,base%m
     swg = e%det_b*base%wg(l) * f(l)
     do i=1,base%pk
       ! cck volume contributions are precomputed
       fk(i) = fk(i) +   &
         swg * base%p(i,l)
     enddo
   enddo
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Now we can exploit the symmetry of the system
   bbk = -transpose(cck)
   !--------------------------------------------------------------------

 end subroutine ldg_locmat_e
 
!-----------------------------------------------------------------------

 !> Compute the local matrices and right hand side (side
 !! contributions).
 !!
 !! The local matrices use a local element based order: the degrees of
 !! freedom are ordered taking first those if <tt>s\%ie(1)</tt> and
 !! then those of <tt>s\%ie(2)</tt>. This function can be used only
 !! for internal sides, so that the second element is always well
 !! defined. For boundary sides, use the function \c ldg_locmat_bs.
 !!
 !! The local matrices correspond to the following terms (notice that
 !! \f$A^s\f$ is not computed at all):
 !! \f{displaymath}{
 !!  A^s_{ij} = 0, \qquad
 !!  B^s_{ij} = \frac{1}{2}\int_s \varphi_j
 !!  \left(n_i\cdot\omega_i\right)\,d\sigma = -C^s_{ji},
 !! \f}
 !! \f{displaymath}{
 !!  C^s_{ij} = -\frac{1}{2}\int_{s} 
 !!  \left(n_j\cdot\omega_j\right)\varphi_i\, d\sigma,
 !!   \qquad
 !!  D^s_{ij} = \int_s {\rm sign}(i,j)\,C_{11}\,
 !!  \varphi_j\varphi_i\,d\sigma.
 !! \f}
 !! Here we use
 !! \f{displaymath}{
 !!  C_{11} = \frac{1}{h}
 !! \f}
 !! and we let (formally)
 !! \f{displaymath}{
 !!  {\rm sign}(i,j) = n_i\cdot n_j.
 !! \f}
 !! The notation \f$n_i\f$ indicates that the normal vector must be
 !! taken as the outward normal of the element where the basis
 !! function with index \f$i\f$ is supported. To ensure consistency of
 !! the formulation even for nonplanar grids, it is essential that we
 !! don't add nor multiply vectors defined on different elements; this
 !! can be obtained by suitably exploiting the (formal) identity
 !! \f$n_{K}=-n_{K^\prime}\f$.
 !!
 !! \bug Most of the computations could be precomputed on the
 !! reference element, however 
 !<
 pure subroutine ldg_locmat_s(      bbs , &
                              ccs , dds , &
                              base,s)
  real(wp), intent(out) :: &
               bbs(:,:) , &
    ccs(:,:) , dds(:,:)
  type(t_base), intent(in) :: base
  type(t_s), intent(in) :: s

  real(wp), parameter :: &
    sij(2,2) = reshape( (/+1.0_wp,-1.0_wp,-1.0_wp,+1.0_wp/) , (/2,2/) )
  integer :: l, ie, i, ii, je, j, jj 
  real(wp) :: h, swg1, swg2,                   &
    no(base%mk,base%ms,2), pb(base%pk,base%ms,2)
  character(len=*), parameter :: &
    this_sub_name = 'ldg_locmat_s'

   ! problem coefficients
   h = 0.5_wp*( diameter(s%e(1)%p) + diameter(s%e(2)%p) )

   ! map the quad nodes on the element quad nodes
   do ie=1,2
     no(:,:,ie) = base%no( : ,                            &
                   base%stab(s%e(ie)%p%pi(s%isl(ie)),:) , &
                   s%isl(ie) )
     pb(:,:,ie) = base%pb( : ,                            &
                   base%stab(s%e(ie)%p%pi(s%isl(ie)),:) , &
                   s%isl(ie) )
   enddo

   !--------------------------------------------------------------------
   ! Second equation: compute ccs and dds
   ccs = 0.0_wp
   dds = 0.0_wp
   do l=1,base%ms
     ! scaled quad weight
     swg2 = s%a/base%me%voldm1 * base%wgs(l) / h
     do ie=1,2
       do i=1,base%pk
         ii = (ie-1)*base%pk + i
         do je=1,2
           ! scaled quad weight
           swg1 = 0.5_wp * base%wgb(l,s%isl(je))
           do j=1,base%mk
             jj = (je-1)*base%mk + j
             ccs(ii,jj) = ccs(ii,jj) -    &
               swg1 * no(j,l,je)*pb(i,l,ie)
           enddo
         enddo
         do je=1,2
           do j=1,base%pk
             jj = (je-1)*base%pk + j
             dds(ii,jj) = dds(ii,jj) +               &
               sij(ie,je)*swg2 * pb(j,l,je)*pb(i,l,ie)
           enddo
         enddo
       enddo
     enddo
   enddo
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Now we can exploit the symmetry of the system
   bbs = -transpose(ccs)
   !--------------------------------------------------------------------

 end subroutine ldg_locmat_s

!-----------------------------------------------------------------------

 !> Compute the local matrices and right hand side (boundary side
 !! contributions).
 !!
 !! This function resembles \c ldg_locmat_s, except that the second
 !! element is absent and in its place the boundary conditions are
 !! used.
 !!
 !! \bug Most of the computations could be precomputed on the
 !! reference element, however 
 !!
 !! \par
 !!
 !! \bug The robin term is still missing.
 !<
 pure subroutine ldg_locmat_bs(      bbs , f1s , &
                               ccs , dds , f2s , &
                               base,s,bs,err)
  real(wp), intent(out) :: &
               bbs(:,:) , f1s(:) , &
    ccs(:,:) , dds(:,:) , f2s(:)
  type(t_base), intent(in) :: base
  type(t_s), intent(in) :: s
  type(t_b_s), intent(in) :: bs
  type(t_bcs_error), intent(out) :: err

  integer :: l, i, j
  real(wp) :: h, swg1, swg2, no(base%mk,base%ms), pb(base%pk,base%ms), &
    xgs(s%m,base%ms), u_b(base%ms), g_b(base%ms), h_b(base%ms)
  character(len=*), parameter :: &
    this_sub_name = 'ldg_locmat_bs'

   ! problem coefficients
   h = diameter(s%e(1)%p)

   ! map the quad nodes on the element quad nodes
   no(:,:) = base%no( : ,                           &
               base%stab(s%e(1)%p%pi(s%isl(1)),:) , &
               s%isl(1) )
   pb(:,:) = base%pb( : ,                           &
               base%stab(s%e(1)%p%pi(s%isl(1)),:) , &
               s%isl(1) )

   !--------------------------------------------------------------------
   ! Boundary data
   xgs = affmap(s,base%xigs)
   bbs = 0.0_wp
   f1s = 0.0_wp
   ccs = 0.0_wp
   dds = 0.0_wp
   f2s = 0.0_wp
   btype: select case(bs%bc)
    case(b_dir)
     u_b = coeff_dir(xgs,-s%ie(2))

     do l=1,base%ms

       ! first equation
       swg1 = 0.5_wp * base%wgb(l,s%isl(1))
       do i=1,base%mk
         do j=1,base%pk
           bbs(i,j) = bbs(i,j) +    &
               swg1 * pb(j,l)*no(i,l)
         enddo
         ! sign change because it goes on the rhs
         f1s(i) = f1s(i) -          &
               swg1 *  u_b(l)*no(i,l)
       enddo

       ! second equation
       ! ccs is zero, since we assume that q is single valued
       swg2 = s%a/base%me%voldm1 * base%wgs(l) / h
       do i=1,base%pk
         do j=1,base%pk
           dds(i,j) = dds(i,j) +  &
             swg2 * pb(j,l)*pb(i,l)
         enddo
         f2s(i) = f2s(i) +        &
             swg2 *  u_b(l)*pb(i,l)
       enddo

     enddo

    case(b_neu)
     ! add here the Robin term
     !g_b = coeff_rob(xgs,s%ie(2))
     h_b = coeff_neu(xgs,-s%ie(2))

     do l=1,base%ms

       ! first equation
       swg1 = 1.0_wp * base%wgb(l,s%isl(1))
       do i=1,base%mk
         do j=1,base%pk
           bbs(i,j) = bbs(i,j) +    &
               swg1 * pb(j,l)*no(i,l)
         enddo
       enddo

       ! second equation
       ! dds is zero, since we assume that u is single valued
       swg2 = 0.5_wp * base%wgb(l,s%isl(1))
       do i=1,base%pk
         do j=1,base%mk
           ccs(i,j) = ccs(i,j) -    &
               swg2 * no(j,l)*pb(i,l)
         enddo
         f2s(i) = f2s(i) +          &
               swg2 *  h_b(l)*pb(i,l)
       enddo

     enddo

    case default
     err%lerr = .true.
     write(err%message,'(A,I5,A,I5,A)') &
       'Unknown boundary condition "',bs%bc,'" on side ',bs%s%i,'.'
   end select btype

 end subroutine ldg_locmat_bs
 
!-----------------------------------------------------------------------

end module mod_ldg_locmat

